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Abstract. 

This  paper  describes  a  solution  to  the  flexural  analysis 
of  a  beam  by  digital  computer.   Subroutines  were  written 
which  calculate  shear,  moment,  slope  and  deflection  at  any 
specified  point  on  a  beam  if  the  loading  is  known.   Shear 
deflection  may  be  included  for  many  support  conditions. 
Neither  elastic  supports  nor  column  effects  are  considered. 
The  subroutines  solve  problems  of  variable  beam  cross- section 
and  loading  by  utilizing  an  extension  of  "McCaulay's  Method", 
a  generalized  step  function. 

The  main  program  provides  input-output  facilities  and 
also  solves  for  Indeterminate  reaotions  on  the  beam. 

Problems  and  their  solutions,  showing  the  versatility 
of  this  program,  are  also  included. 
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1.   Introduction. 

The  routine  solution  to  engineering  problems  by  digital 
computer  appears  to  be  limited  only  by  the  scope  of  the 
available  programs.  A  flexible  program  permitting  the  solu- 
tion of  routine  problems  in  flexural  analysis  of  beams  would 
seemingly  fill  a  need  in  the  field  of  digital  computer  ap- 
plications to  Mechanics.   The  present  work  is  devoted  to  such 
an  idea. 

Five  subroutines,  "LOAD",  "SHEAR",  "MOMENT",  "SLOPE", 

and  "DEFLECT",  were  written  which  will  calculate  the  load, 

shear,  moment,  slope,  and  deflection,  respectively,  at  any 

indicated  point  on  a  beam  if  the  indeterminate  quantities 

are  calculated  elsewhere.   The  use  of  generalized  step  functions 

in  this  project  allows  for  piecewise  linear  variation  of 

distributed  loading,  bending  compliance,  — ,  and  shear  com- 
ic 
pliance,  -jg. 

The  main  program,  "BEAM3",  first  generates  and  then 
solves  a  set  of  simultaneous  linear  equations,  calculating 
the  indeterminate  quantities  of  the  beam.   The  above  mention- 
ed subroutines  are  then  utilized  to  calculate  the  remaining 
flexural  quantities.   This  program  also  provides  input-output 
facility. 


2.   Definition  of  the  problem. 

The  problem  of  structural  beam  analysis  is  one  in 
which  the  external  loading  applied  to  the  beam  is  well  de- 
fined, as  is  the  geometric  configuration  and  the  nature  of 
the  beam  itself.   Where  information  about  such  a  beam  is  re- 
quired, this  information  will  most  likely  fall  into  the 
following  categories: 

1.  Internal  shear  forces,  V,  and  moments,  M. 

2.  Slope, £. 

3.  Deflection,  y. 

k.      Redundant  reactions,  either  forces,  Ry,  or  moments, 
RMZ. 

"Flexural  analysis"  is  used  here  as  meaning  the  process  by 
which  this  information  is  generated  from  the  known  data. 

In  attacking  this  problem  of  flexural  analysis,  the 
intention  was  to  limit  the  scope  as  little  as  possible. 
However,  several  customary  simplifications  and  limitations 
were  made.   First,  elastic  supports  and  beam  column  effects 
were  not  considered.   Limitations  concerning  the  type  and 
number  of  external  loads  and  the  geometry  of  the  beam  will 
be  discussed  later.   Also  the  basic  relations  which  follow, 
inherently  restrict  the  deflection  of  the  beam  to  small  values 
and  the  internal  stresses  to  values  below  the  elastic  limit 
of  the  structural  material.   Shear  deflection,  which  is 
generally  neglected  in  a  problem  of  this  type,  has  been 
included  in  some  problems. 
a.   Bending  deflection  (y1) . 

The  Bernoulli-Euler-Navier  equation  governing  the  elastic 
curve  of  a  beam  due  to  bending,  consistent  with  the  sign 
convention  shown  in  Appendix  B  is: 


—  ■  curvature  «  (d  y/dx  )  [l^Cdy/dx)  J 

2    2 

and  we  use  the  approximation,  curvature  £=z  d  y/dx  . 

The  following  equations,  derived  from  the  above  equation 
and  shown  in  Timoshenko  /5/,   with  differences  due  to  the  chosen 
sign  convention,  describs  the  flexure  of  a  beam  (See  Appendix  B 
for  notation) . 

W  =  q(x)   , 

V  »  /wdx   , 

M  «  yVdx   , 

*-£  -A*  • 

b.   Shear  deflection  (y2). 

The  equation  governing  the  shear  deflection  of  a  beam, 

also  given  in  Timoshenko  /5/»  compatible  with  the  stated 

sign  convention  becomes: 

dy2  m         K_V 
dx        AG 
It  is  evident  that: 

?2      =j\   AG/  dx 

o.   Total  deflection  (y). 

The  equations  used  in  this  thesis  are  the  result  of  the 
sum  of  the  above  shear  and  bending  deflections.   They  are: 

W   -  q(x)   , 


Numbers  so  indicated  refer  to  the  bibliography, 


V 

=     [    wax    +  EFj* 

M 

=       P  Vdx     ♦    ZMj* 

e 

fx  M                      KV 
"      \      EI  dX     -      AO 

-JO 

Y 

■     1    0  dx    +    c2** 

Jo 

♦  c  ## 


Note  that  step-discontinuities  may  appear  in  all  of  the 
above  equations,  except  in  the  equation  for  deflection.   These 
discontinuities  limit  the  methods  that  can  be  used  in  handling 
the  equations. 

Three  methods  that  could  be  used  to  solve  the  equations 
are  a  pure  analytical  approach,  a  step-by-step  numerical 
integration  scheme,  and  the  "curly  bracket"  method.   Of 
these  methods,  the  analytical  approach  is  generally  restrict- 
ed to  a  specific  problem  and  is  unsuited  for  adaptation  to 
digital  computer  programming.   The  numerical  integration  scheme 
can  be  used  in  the  computer  solution  of  flexural  analysis 
and  may  have  been  programmed;  however,  our  search  of  literature 
has  not  uncovered  any  such  program.   The  "curly  bracket" 
method  was  chosen  for  use  in  this  thesis  because  it  was  felt 
that  it  would  result  in  a  more  exact  solution  requiring  less 

♦In  these  equations,  ZTF1  and£M   represent  the  sum  of 
the  point  forces  and  moments,  respectively,  applied  to  the 
beam. 

♦•C^  and  C2  are  constants  of  Integration. 


computer  time,  particularly  if  results  were  required  at  only 
a  few  points. 


3.   "Curly  bracket"  method 

The  "curly  bracket"  method  is  based  upon  the  properties 
of  a  convenient  notation  defined  as  follows: 

{x-a}         =0  if  x<a 

*        (x-a)  if  x>a 

where  n  is  an  integer^O. 

Appendix  A  contains  further  properties  of  expressions  contain- 
ing curly  brackets.   The  notation  is  generally  credited  to 
McCaulay  /l/,  but  similar,  though  less  useful  notations  can 
be  traced  to  A,F6*ppl  and  St.  Venant.   Also,  the  ideas  in- 
volved are  related  to  the  use  of  the  Laplace  Transformation 
with  the  result  that  this  method  or  approximations  to  it  are 
reinvented  frequently  /8/  /9/. 

Using  this  method  allows  one  equation,  valid  over  the 
entire  length  of  the  beam,  to  be  expressed  for  each  of  the 
basic  equations  of  flexure.   The  particular  value  of  this 
method  is  that  it  permits  the  integration  of  these  equations 
to  be  carried  out  easily  and  has  the  virtue  of  adjusting  the 
constants  of  integration  most  conveniently. 

In  addition,  for  use  with  a  digital  comouter,  the  "switch- 
ing" property  of  the  curly  bracket  can  be  easily  handled  in 
FORTRAN  by  use  of  an  "IF"  statement. 

A  simple  example  should  convince  most  readers  of  the 
convenience  of  this  method.   Examine  the  following  problem, 
which  consists  of  a  uniform  beam,  simply  supported  at  each 
end  with  a  concentrated  load  at  the  center  of  the  soan. 


a 


' 


2  2 

Fig.   1 
The  following  equations  lead  to  the  "solution"  of  this 
problem. 


V   - 
M   - 

EI 
Ely 


2    *  * 


;■ 


*2    V  ' 


12    "   g 


Cx  +  /c-r;- 


The  only  unknown  value,  C-   in  the  equations  can  easily  be 
calculated,  using  the  boi     '  condition,  y=0  at  x=2a. 


0-, 


-8' 


a. 


. 


gsuiting  equation  defining  deflection 
along  the  beam  is: 


t  any  point 


y  = 


• 


p 


3        -„a^3 
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As  McOaulay  used  thi     Shod,  it  i :  suitable  only  for 

beams  of  uniform  crc    section.   Quinlan  /3/  -     sted  a  way 
of  excan  5    the  method  jo  include  beams  with  piecewise  con- 


stant  variation  of  bending  compliance,  •==.   In  this  thesis 
the  method  has  been  further  expanded  to  include  beams  with 
piecewise  linear  variation  of  both  bending  compliance  and 
shear  compliance,  -— . 

All 
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k.      Implementation. 

Although  the  formulas  shown  in  Appendix  A  could  be  ex- 
panded further,  a  piecewise  linear  variation  of  distributed 
loads,  bending  compliance,  and  shear  compliance  was  consider- 
ed capable  of  providing  sufficient  accuracy  and  flexibility 
for  the  resulting  program. 

Fig.  2(a)  illustrates  a  distributed  load  on  a  beam 
whose  origin  is  at  x=0.   Fig.  2(b)  shows  an  element  of  the 
distributed  load. 


Fig.    2(a) 


T 

Qa 


A — * 
—  B 


-*» 


% 


Fig.    2(b) 
The  analytical  expression  for   this   element  is: 

Q(x)      =      Qajx-Aj°      *      V^x-A}1-      %  {x-BJ  °  -  %!**    {x-B> X 

B— A  B— A 


The  same  relationship  is  used  for  an  element,  f(x),  of  bend- 
ing compliance  and  an  element,  g(x),  of  shear  compliance  upon 
substituting  f  for  Q  and  g  for  Q  respectively. 

The  basic  equations  of  flexure  can  now  be  expressed  in 
the  same  manner  as  they  were  utilized  in  programming.   They 


are; 


NOL 

W(x)     -     2    <M») 

i=l     x 


(1) 


/JC  NOF  n  NOR 

V(x)      «  W(x)dx     t      S     Fy    fx-xf  r      ♦      E»,    fr-*r  1 


(2) 


where  F   are  upward,  known  forces  at  x=xf  and  R   are  up- 

yl  1      yi 

ward,  unknown  reactions  at  x=x 


NOM 


NORM 


M(X)        =     (       V(X)dX       +        2      M*      fX~Xm    )°       ♦         S        RM        JX-X  V 

(3) 

where  M   are  clockwise,  known  moments  at  x=x   and  RM_  are 


mi 


clockwise,  unknown  reactive  moments  at  x^x 


rm< 
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dy 
dx 


NOI 
M(x)  2  f.(x) 
i=l 


dx  - 


NOK 
V(x)  Z  g4(x) 
i=l   X 


♦   C 


W 


where  C,  is  the  constant  of  integration. 

Appendix  A  indicates  how  the  product  and  the  integration  of 
the  product  of  the  curly  bracket  symbols  are  formed. 


rxrxr   a  NOI       -i        rx-    NOK    n 


dx  ♦  C,x  ♦  C2 


(5) 


where  C^  is  the  constant  of  integration. 

Equations  (1)  through  (5)  when  properly  coded  in  FORTRAN 

2The  numbers,  NOL,  NOF,  NOR,  NORM,  NOI  and  NOK,  appear- 
ing above  the  summation  signs  are  program  input  parameters 
indicating  the  number  of  terms  for  each  problem. 
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form  the  core  of  the  five  basic  subroutines i      LOAD,  SHEAR 
MOMENT,  SLOPE  and  DEFLECT.   It  must  be  noted  before  discuss- 
ing the  subroutines  that  equations  (2)  through  (5)  involve 
the  unknown  values  0^,  02,  Ry  »  and  RM   ,  where  there  may  be 
as  many  as  ten  values  each  of  Ry  and  RMz„   The  calculation  of 
these  values  is  done  in  the  main  program  BEAM3  and  will  be 
discussed  in  section  6. 
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5.   The  basic  subroutines. 

Subroutines  LOAD,  SHEAR  and  MOMENT. 

The  above  named  subroutines,  as  their  names  imply,  are 
employed  to  calculate  the  total  distributed  load,  shear  and 
moment,  respectively,  at  any  desired  point  on  the  beam.   The 
expanded  forms  of  equations  (1),  (2)  and  (3)  were  used  in 
writing  these  subroutines.  Utilizing  the  distance,  x*,  as 
the  entering  argument,  the  subroutines  select  and  sum  the 
appropriate  terms  in  the  corresponding  equation  by  inspect- 
ing all  loading  terms  and  discarding  those  in  which  the  value 
within  the  curly  bracket  is  negative.  Actually  the  load, 
shear  and  moment  are  calculated  a  small  increment,  £   ■   10   , 
to  the  right  of  x^  in  order  to  obtain  the  right  side  value 
of  any  step-discontinuities  oocuring  at  the  point,  x,  (See 
flow  charts  on  pages  13,  1^  and  15). 

Subroutine  SLOPE. 

This  subroutine  is  equation  (4-),  excluding  the  constant 
of  integration,  C^  in  FORTRAN  language  (See  flow  chart  on 
page  17).   Note  that  the  equation  can  be  expressed  as  a 
summation  of  terms,  which  are  constants  times  the  multiplica- 
tion of  two  step-functions  of  the  form  C  £x-a] n [x-b} m,  or  the 
summation  of  the  integral  of  the  same  type  of  terms,  where  n 
equals  either  0  or  1. 

Function  subroutines,  FTON  and  FT1N  (see  Appendix  A), 
will  calculate  the  value  of  £x-aj  [x-bjm  and  [x-aj  {x-b}  n, 
respectively.   Function  subroutines,  ENTON  and  ENT1N  (See 
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Flow  Chart   of   LOAD 


SUBROUTINE 
LOAD(X.W) 


o 

EPSX=1Q~ 

•< 

' 

¥=0 

-=<X+KP3X-B(ir>^ 


0 


r—V  X + EP  8X- A  ( ij^ 


Note:   Numbers  preceded  by  a  number  sign,  #,  refer  to 
the  statement  number  in  the  listing  of  LOAD. 
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Flow  Chart  of  SHEAR 


SUBROUTINE 
SHEAR(X.V) 

■ 

i 

▼-0 

J 

I 

EPSX-10"8 

N-l 


VV-FY(N)    | 


VV«  #4 


vv-  #5 


v-v+w 


1  I-."11 


0 


M-l 


+ 

'    X*EPSX 
\-XR(M)     . 

V- 

1 

■ 

VV-RY(M) 

o 

■ 

r 

v-v*vv 

' 

r 

M-H*l 

■ 

' 

Mote:      Numbers   Dreceded  by  a  number   sign,    #,    refer   to 
the   stufpnient   number   In    tine    listing  of   SHEAR. 
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Flow  Chart  of  MOMENT 


SUBROUTINE 
MOMENT   (X.AM) 

1 

' 

EPSX-10-8 

AM-0 

RETURN 


: 


Not*:   Numb?rs  preoe4e£  by  a  number  sign,  #,  refer  t 
statament  number  In  the  listing  of  H  . 
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Appendix  A),  will  calculate  the  first  integral  of  these 
functions.   Equation  (*0  can  now  be  expressed  as  the  summation 
of  constants  multiplied  by  the  appropriate  function  subroutines. 

SLOPE  calculates  all  combinations  of  M(x)  and  f(x)  by 
stepping  through  all  the  f(x)  terms  and  varying  all  the  terms 
in  the  moment  equation  at  each  step.   It  then  steps  through 
all  the  g(x)  terms  varying  all  the  terms  in  the  equation  for 
shear  at  each  step.   The  sum  of  these  terms  is  the  value  of 
the  slope,  minus  the  integration  constant,  C^,  at  the  entered 
distance,  x. 

Subroutine  DEFLECT. 

At  this  point  the  author  wishes  to  point  out  the  similar- 
ity of  the  equations  for  slope  and  deflection.   Since  the 
integration  constants,  Cn  and  Cg,  are  again  excluded  from  the 
subroutine,  and  the  function  subroutines,  DITON  and  DIT1N 
(See  Appendix  A),  calculate  the  second  integral  of  the  mul- 
tiple step-functions,  subroutine  DEFLECT  is  identical  to 
SLOPE  except  for  substituting  ENTON  for  FTON,  ENT1N  for  FT1N, 
DITON  for  ENTON  and  DIT1N  for  ENT1N.   The  result  is  a  sub- 
routine that  will  oalculate  the  deflection,  minus  the  con- 
stants of  integration,  at  any  required  point  on  the  beam. 
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6.   Program  BEAM3. 

The  main  program,  BEAM3,  oan  be  divided  into  four  dis- 
tinct steps.   These  steps  are:   reading  the  input  data,  cal- 
culation of  indeterminate  quantities,  calculation  of  the 
flexural  quantities,  and  the  output  of  the  solution. 

The  program  first  reads  the  input  data  which  includes 
the  length  of  the  beam,  all  external  loading,  the  location  of 
all  reactions,  the  flexural  properties  of  the  beam,  and  a  group 
of  parameters  giving  the  number  of  each  type  of  load  or  re- 
action.  It  also  provides  for  reading  a  predetermined  value 
of  deflection,  PDEFT,  at  each  reaction  and  a  predetermined 
value  of  slope,  PSLP,  at  each  moment  reaction. 

The  program  then  proceeds  to  solve  for  the  indeterminate 
quantities  of  the  beam  by  utilizing  the  two  equations  of 
static  equilibrium  plus  equations  obtained  from  the  known 
boundary  condition  at  each  reaction  or  reactive  moment.   The 
first  two  equa;ions  are  simply  the  summation  of  the  forces 
on  the  beam  equals  zero  and  the  summation  of  the  moments 
about  the  right  end  of  the  beam  «quals  zero.   The  remainder 
of  the  equations  are  obtained  from  the  fact  that  the  slope 
at  each  reactive  moment  or  the  deflection  at  each  reaction 
is  known.   Using  this  fact  in  the  appropriate  one  of  either 
equation  (*0  or  (5)  results  in  a  set  of  linear  equations  for 
the  indeterminate  quantities. 

The  equations  thus  obtained  are  generated  in  a  matrix 
form  that  is  compatible  with  subroutine  GAUSS2  (See  Appendix  E). 
This  subroutine  is  then  called  to  solve  the  equations  simulta- 
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neou8ly.   If  the  matrix  Is  singular,  GAUSS2  exits  to  the  call- 
ing program  which  in  turn  prints  "MATRIX  SINGULAR"  and  then 
stops.   If  not,  the  subroutine  solves  the  equations  and  re- 
turns the  calculated  values  to  the  main  program.   The  most 
serious  limitations  of  this  program  are  Imposed  by  the  gen- 
eration or  solution  of  the  simultaneous  equations,  therefore 
it  is  appropriate  to  include  and  explain  several  of  them  at 
this  point. 

First,  at  least  one  reaction,  even  if  it  is  later  forced 
to  equal  zero,  must  be  included  with  the  input,  or  a  singular 
matrix  will  be  generated. 

Including  shear  deflection  in  the  solution  must  also  be 
greatly  limited  at  this  point.   Since  in  general  the  shear 
includes  step  discontinuities,  they  will  also  appear  in  the 
slope  if  shear  deflection  is  included.   The  program  generates 
an  equation  for  slope  at  each  moment  reaction.   If  by  chance, 
either  a  reaction  or  point  force  occurs  at  the  same  point  as 
a  moment  reaction,  the  slope  at  this  point  is  double  valued. 
The  program  is  unable  to  distinguish  which  of  the  two  values 
is  required  except  when  it  occurs  at  the  extreme  ends  of  the 
beam.   For  example  a  beam  "built  in"  at  both  ends  is  accept- 
able.  If  the  beam  is  "built  in"  at  only  one  end,  it  must  be 
solved  in  a  manner  such  that  the  "built  in"  is  at  the  origin 
of  the  beam. 

The  program  then  divides  the  beam  into  a  number  of 
increments,  NOP-1,*  and  calls  the  subroutines,  LOAD,  SHEAR, 

*NOP  is  another  parameter  read  in  at  the  start  of  BEAM3. 
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MOMENT,  SLOPE  and  DEFLECT  to  calculate  the  flexural  quantities 
at  each  of  these  increments. 

The  output  consists  of  printing  the  distance  from  the 
left  end  of  the  beam  to  the  point  in  question  and  the  values 
of  the  load,  shear,  moment,  slope  and  deflection  at  that 
point.   The  values  of  all  of  the  indeterminate  quantities  are 
also  printed  out. 
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Fig,  7      Simplified  Flow  Chart  of  BEAM3 


START 


READ   INPUliP) 


CALCULATE  TWO 
ROWS  OF  MATRIX, AA, 
SATISFYING   STATIC 
EQUILIBRIUM 


\f 


CALCULATE  REMAINDER 
OF  MATRIX, AA, 
SATISFYING  ELASTIC 
EQUILIBRIUM 


i_ 


CALL 
GAUSS2 


CALCULATE 
W,V,M.e,Y 
AT    "NOP5'   POINTS 


MATCIX  |     STop 


SINGULAR 
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7.   Testing  of  Program. 

The  testing  of  the  program  and  the  subroutines  could, 
in  theory,  continue  indefinitely  because  of  the  unlimited 
number  of  beam  configurations.   The  author  has  attempted  to 
show  (See  Appendix  D)  solutions  to  test  problems  of  simple 
configurations  and  constant  bending  and  shear  compliance  to 
compare  with  classical  solutions.   Then  a  statically  deter- 
minate problem  with  a  varying  bending  compliance  was  solved 
and  compared  to  a  solution  obtained  from  a  numerical  in- 
tegration scheme.   The  most  rigid  test  for  the  program  was 
in  constructing  a  problem  for  a  beam  of  unlikely  configuration 
(See  Problems  6  and  7  of  Appendix  D),  and  solving  the  flexural 
analysis  of  this  beam  from  both  ends.   Note  that  the  only 
significant  difference  between  the  two  solutions  is  the  opposite 
signs  for  shear  and  slope,  which  would  of  course,  be  expected. 
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8.   Conclusions  and  Recommendations. 

The  thesis  shows  a  method  of  including  pieoewlse  linear 
variation  of  both  bending  compliance  and  shear  compliance  in 
the  flexural  analysis  of  beams  by  extending  "McCaulay's  Method". 
Adapting  this  idea  to  the  digital  computer  resulted  in  a 
program  capable  of  analyzing  a  very  general  type  of  beam. 
However,  had  time  permitted,  several  changes  would  have  been 
made  to  the  program.   First,  internal  monitoring  of  the  input 
data  would  be  an  asset  to  any  user.   Also,  using  the  same 
general  method  of  attack,  this  program  would  have  been  ex- 
panded to  include  bending  compliance  and  shear  compliance  of 
a  piecewise  parabolic  nature. 

Other  additions  to  the  program,  much  larger  in  scope, 
would  be  including  elastic  supports  and/or  beam  column  effect. 

It  would  be  interesting  to  compare  the  results  of  this 
program  to  those  obtained  from  a  program  utilizing  the  numer- 
ical integration  method.  Also,  a  more  exhaustive  study  of  the 
program  could  be  made  to  determine  accuracy  of,  and  time 
requirement  for,  solutions. 
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APPENDIX  A 
"Curly  bracket"  method,  definitions  and  formulas 
We  define: 

1.   {x-a}n    =     (x-a)n    ,    x>a 
=     0       ,    x<a 
where  n  is  an  integer ^0 
From  this  definition,  it  is  possible  to  establish  the 
following  formulas .   Although  it  would  be  possible  to  express 
these  results  in  terms  of  curly  brackets,  in  most  cases  the 
resulting  expressions  would  be  rather  complicated.   The  forms 
exhibited  below  have  the  advantage  of  being  most  readily 
conformable  to  the  decision  commands  available  in  FORTRAN. 


f  [x-a]". 
Jo 


3.       [x-a}°{x-b3n 


*.       $x-a}1£x-b]n 


5.   JX)x-a}0[x-b}ndx 


(x-a) 
n+1                » 

x>  a 

0 

Otherwise 

(x-b)n 

x>  a  and  x>b 

0 

Otherwise 

(x-b)          + 

(b-a)(x-b)n 
and 

,    x>a 

x>b 

0 

Otherwise 

t      mn+1 
(x-b) 

n+1 

(a-b)n+1 
n+1 

x>a>b 

(x-b) 
n+1 

» 

x>b>a 

,      Otherwise 
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6.    I    [x-aj     [x-b]    dx       =       -^ 
Jo 


(x-b)  +  (D-a)(x-b) 


n+1 

(a-b) 
*      (n*l)(n+2)  '      x>a>b 

.      U^^2   ,      (D-aHx-D)"*1        x>b>a 

n+2  n*l 


,    Otherwise 


7-  C£  B-a3°^3n^ 


x-b )  (x-a) (a-b) 

n*l)(n+2)  (n+1) 


x>a>b 


f        v.xn  +  2 

(a-b) 
(n+l)(n+2) 

n+2 
(x-b) 

(n+l)(n+2)  ,  x>b^a 

0  ,  Otherwise 


rX     /X 


8-Ii  $-?Mn**' 


/-       ^xn*3  /v       w      ^n*2 

(x-b)  (b-a) (x-b) 

(n+2)(n+3)  (n+l)(n+2) 
*-a, (a-b)  2(a-b) 


*      (A+3);n+2)  (n+l)(n+2)(n+3)    ,      x>a>b 

(x-b)  (b-a) (x-b)  v^ 

(n+2)(n+3)  (n+l)(n+2) 

«  0  ,    Otherwise 

Function   Subroutines   FTON,    FT1N,    ENTON,    ENT1N,    DITON,    DIT1N 

The   above  named  function   subroutines  were  written  and 
utilized  to  evaluate   each  of  the  preceeding   formulas  numbered 
3  through  8.      The  values  x9    a,    b     and  n  were   used  as  entering 
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arguments   to   these   functions.      The   function   then  became   the 
variable  as   follows: 


FTON  =  £x-a}°£x-b}n 

FT1N  -  jx-a}     £x-b] 

EN  TON        =        (      £x-a}  °  [x-b]  ndx 
^o 

EN  TIN       =        V      jx-a]}    [x-b}ndx 
Jo 

DITON        -        Tl       [5~aj°  fe-b?na^  dx 
DIT1N        =        (*(*  [^-a] 1  [^-bjnd^  dx 


As  examples,    flow  charts   for  FTON  and  ENTON   are   shown  in 
Figs.    8  and  9, 
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Fig.  8      Flow  Chart  of  FTON 


function 
fton(x,a,b,n) 


FTON»0 


FT0N=(X-B) 


N 


RETURN 
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Fig.  9       Flow  Chart  of  ENTON 


FUNCTION 
ENTON    (X,A,B,N) 


FN=N 


EN  TON =0 


(X  B)N*1    (A  R)N  +  1 
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APPENDIX  B 
Sign  Convention  and  Notation 
Sign  Convention 

The  following  sign  convention  has  been  employed: 

y 

x  y  ■  deflection 


a 


f    X 


©  ■   slope  ■  _Z 

dx 


M 


M 


bending  moment 


vtC3 


V     ■     shear  force 


W 


rm 


Note 


W 


distributed  loads 


:   1.   Point  forces,  Fy,  and  point  reactions,  R  ,  are 
positive  in  the  positive  y  direction. 

2.   Point  moments,  M2,  and  point  reactive  moments, 
RM2,  are  positive  in  a  olockwise  sense. 


Notation; 

THESIS 

PROGRAM 

C 

Length  of  beam 

W 

w 

Total  distributed  load  at  a  point 

V 

V 

Shear 

M 

AM 

Moment 

e 

Yl 

Slope 

y 

y 

Deflection 
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*y 

FY 

Ry 

RY 

Mz 

AMZ 

RMZ 

RAMZ 

Q(x) 

Qa 

AQ 

% 

03 

f(x) 

EA 

EB 

Point  load  or  foroe 

Point  reaction 

Point  moment 

Point  reactive  moment 

Element  of  distributed  load 

Distributed  load  (left  end) 

Distributed  load  (right  end) 

Element  of  bending  compliance 

Bending  compliance      ax 
(left  end  of  element) 

Bending  compliance 
fright  end  of  element) 
g(x)  Element  of  shear  compliance 

NAG; 
Shear  compliance 
(left  end  of  element) 
Shear  compliance 
(right  end  of  element) 
Distance  to  Fy 

Distance  to  Ry 

Distance  to  Mz 

Distance  to  RMZ 

Distance  to  Ga 

Distance  to  G^ 

Distance  to  Ea 

Distance  to  E^ 

Distance  to  Q_ 

Distance  to  Q,^ 

Number  of  distributed  loads  (25)  ^ 

3 
Numbers  following  explanation  of  terms  NOL — NOK  indicate 

the  arbitarily  chosen  maximum  value  of  each  for  this  program. 
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GA 

GB 

xf 

XFO 

*r 

XR 

*m 

XM 

Xm 

RXM 

Ag 

AG 

Bg 

BG 

Ae 

AE 

Be 

BE 

A 

A 

B 

B 

NOL 

NOL 

NOF 

NOF 

NOR 

NOR 

NOM 

NOM 

NORM 

NORM 

NOI 

NOI 

NOK 

NOK 

PDEFT 

PDEFT 

PSLP 

PSLP 

NOP 

NOP 

Number  of  point  forces  (100) 

Number  of  point  reactions  (10) 

Number  of  point  mom3nts  (50) 

Number  of  point  reactive 

moments  (10) 

Number  of  bending  compliance 

elements  (25) 

Number  of  shear  compliance 

elements  (25) 

Predetermined  deflection 

at  a  reaction 

Predetermined  slope  at  a 

reactive  moment 

Number  of  points  at  which 

the  flexural  quantities  will 

be  calculated. 


32 


APPENDIX  0 
General  Information  Concerning  Use  of  Program 

The  composite  program,  BEAM3  plus  subroutines,  was 
written  in  FORTRAN-60,  compiled  and  tested  on  a  Control 
Data  Corporation  1604  digital  computer.   It  consists  of  412 
FORTRAN  statements  requiring  about  580  cards.   The  computer 
storage  requirement  is  about  12,500  cells,  however,  this 
number  can  be  reduced  by  changing  the  dimension  statements 
which  would  of  course,  reduce  either  the  maximum  number  of 
external  loads,  elements  of  shear  compliance  or  bending  compli- 
ance, or  number  of  points  at  which  flexural  quantities  could 
be  calculated. 

The  FORTRAN  statements  used  in  this  program  are;   GO  TO, 
computed  GO  TO,  IF,  STOP,  DO,  CONTINUE,  FORMAT,  READ,  PRINT, 
WRITE  OUTPUT  TAPE,  FUNCTION,  SUBROUTINE,  RETURN,  CALL,  DIMENSION, 
COMMON,  and  END, 

There  are  no  conversion  constants  incorporated  in  this 
program.  As  such  all  input  must  be  in  a  oonsistant  set  of 
units. 

"FORMAT"  for  Input  Data 

The  input  data  for  a  problem  is  read  by  the  program 
with  a  FORMAT  of  14  for  all  "fixed  point M  variables  and  F20.0 
for  all  floating  >oint  variables.   The  arrangement  of  the 
data  on  the  input  cards  is  shown  as  follows $ 
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Group 
no. 

No 
i 

8  of  cards 
n  group 

1 

1 

2 

1 

3 

NOL 

4 

NOF 

5 

NOR 

6 

NOM 

7 

NORM 

8 

NOI 

9 

NOK 

Fields 

F20.0 

814 

4F20.C 

2F20.0 

3F20oO 

2F20.0 

3F20.0 

4F20.0 

4F20 . 0 


Entries 
G 

NOP,  NOL,  NOF,  NOR 
NOM,  NORM,  NOI,  NOK 
A,  B,  Aa,  Ab 

*f,   Fy 

xr,  0.0,  PDEFT 

xm»  Mz 

xm»  0.0,  PSLP 

Ae»  Be»  Ea»  ^ 

Ag»  Bg»  aa»  °b 


As  an  example  the  input  data  cards  for  test  problem  #1 
in  Appendix  D  was  as  follows: 


CARD 

Is 

100 

CARD 

2s 

21   1   0   2  0 

2   1 

0 

CARD 

3s 

0. 

100 

-500. 

-500, 

CARD 

4s 

0. 

0. 

.0 

CARD 

5s 

100 

0. 

oO 

CARD 

6s 

0. 

0. 

0. 

CARD 

7s 

100 

0. 

0. 

OARD 

8s 

0, 

100. 

.000000001 

.000000001 
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APPENDIX  D 
Test  problems  and  solutions 
Test  problem  1 

This  problem  consists  of  a  beam  of  uniform  cross-section, 
"built-in"  at  both  ends,    and  loaded  with  a  uniformly  dis- 
tributed weight  of  500#/in.    (See   Fig.    10) 


n"  10      #Fn.2 


K 

AG 


.   0 


Fig,    10 


Test  problem  2 

This  problem  is  identical  to  problem  1  except  that  the 
slope  at  either  end  of  the  beam  Is  forced  to  values  other 
than  zero. 


9 


left  end 


002 


right  end 


002 
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Test  problem  3 

This  problem  consists  of  a  beam  of  uniform  cross-section, 
"built-in"  at  one  end  and  simply  supported  at  the  other,  and 
loaded  with  a  uniformly  distributed  weight,   (See  Fig.  11) 


V 


s 


c 


4  -  10#/in. 


ii  i  i  r 


s  N  =-t  - 


40  in, 


I 

EI 

L. 

AG 


1332x10 


-7  1 


#in.2 


Fig.  11 

Test  problem  k 

This  problem  is  identical  to  problem  number  3  except 
shear  deflection  is  included  in  the  solution. 
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Test  problem  5 

This  problem  consists  of  a  homogeneous  circular  shaft, 
three  Inches  in  diameter  at  the  center  and  tapered  to  one 
and  a  half  inches  in  diameter  at  both  ends.   The  total  weight 
of  the  shaft  is  100#  and  a  uniformly  distributed  load  total- 
ing 200#  is  placed  on  the  central  section.   The  shaft  is 
supported  by  a  uniformly  distributed  reaction  at  each  end. 
Points  A  and  B  are  points  of  zero  deflection.   This  problem 
was  approximated  for  the  program  as  shown  in  Fig.  12. 

8.33  Vm. 


Fig.  12 


Weight  of  shaft: 


0"  -  8" 
8"  -   24" 
24"  -   36" 


-0.5315#/in. 

-0.5315   to  -2.o81#/ln. 

-2.081#/lno 

Bending   compliance  I— J  j 

0.1420x10" ,  l/#-in?        , 

0.1420x10"°   to  0.0583x10"° 

0.0583x10"?   to  0.0281xl0"9 

0.0281x10"?  to  0.0152x10"° 

0.0152x10*°   to  0.0089x10"° 

0.0089x10'°  l/#-in? 


(linear) 


0"  -  8H 
8"  -  12" 
12"  -  16" 
16"  -  20" 
20"  -  24" 
24"  -   36" 


i/#-mf 

l/#-in?  (linear) 

l/#-ing  (linear) 

l/#-in?  (linear) 


(linear) 
Clin 


Shear  compliance  VAG/    ~  ® 
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Test  problem  6 

This  hyperstatic  problem  consists  of  a  beam  with  a  very 
unlikely  configuration  as  shown  in  Fig.  13.   The  bending  com- 
pliance varies  as  shown  in  Fig.  lk„    Shear  deflection  has  been 

K 
omitted,    i.e.,  -^  ■  0. 
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Fig.    Ik 

Test  problem  7 

This  problem  is  identical  to  problem  number  6  except  the 
beam  was  turned  end  for  end. 
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APPENDIX  E 
Subroutine  GAUSS2 

This  routine  is  a  modified  form  of  the  GAUSS2  used  by- 
Co  B,  Bailey  /7/  in  his  program,  LINEQN,  for  solving  simul- 
taneous linear  equations  by  Gaussian  elimination,,  Bailey 
was  solving  the  matrix  equation,  Ax^B,  with  a  possibility 
of  50  vectors  B  for  each  matrix  A,  but  this  program  requires 
only  one  vector  B  with  each  matrix  A0  Also  since  the  maximum 
number  of  equations  in  BEAM3  is  22,  the  size  of  matrix  A  was 
reduced  from  100x101  to  22x230 

At  each  step  of  the  elimination,  the  value  of  the  diagonal 
element  is  compared  to  a  defined  MzeroM.   If  it  is  smaller 
than  this  quantity,  the  matrix  is  considered  singular  and 
an  error  output  is  returned  to  the  calling  program.   This 
is  the  "matrix  singular"  test  used  by  BEAM30 
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APPENDIX  F 
Listing  of  program  and  subroutines 

Name  Page 

BEAM3  ^8 

LOAD  5^ 

SHEAR  55 

MOMENT  56 

SLOPE  58 

DEFLECT  61 

FTON  and  FT1N  6k 

ENTON  and  EN TIN  65 

DITON  and  DIT1N  66 

GAUSS2  6? 
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PROGRAM    BEAM3 

DIMENSION    W(  1001) ,V( 1001  ) , AM ( 100 1 ) , Y 1 ( 1 00 1 ) , Y { 1 001  ) , XX(  10 01 )  , 
1     A(25),B(25) ,QA(25),QB(25) ,XF0(100),FY( 100) ,XR( 10},RY{ 10) , 
2XM(50) ,AMZ(50),RXM( 10) ,RAMZ( 10 ) , AE( 25 ) , BE (25 ) , E A ( 25 ) , EB { 2 5 ) , 
3AA(22,23),XA(22) ,PDEFT(1  0),PSLP(10) 
h    , ITITLE( 10) ,GA(25),GB(2  5),AG(2  5) ,BG(25) 

COMMON    A,B,QA,QB,XFO,FY, XR, RY , XM, AMZ, RXM , RAMZ, AE, BE, EA,S3,N0L, 
1     NOR, N0F,N0M, NORM, NOI ,N3K, GA,GB, AG ,BG 

READ       70, C 

70  FORMAT    (F20.0) 

READ    71 ,N0P,N0L,N0F,N0R, NOM, NORM, NOI , NOK 

71  F0RMAT(8IU) 
IF(NOL)     51,51,50 

50  READ    72, (A( I  ),B(I ),QA(I)  ,QB( I) ,1=1, NOD 

72  F0RMATUF20.0) 

51  IF(N0F)53,53,52 

52  READ  73, (  XFO{ N) ,FY (N) ,^= 1 , NOF) 

73  FORMAT  (2F20.0) 

53  IF(N0R)55,55,5U 

54  READ  74, (XR(M),RY(M) ,PDE FT ( M ) , M= 1 ,N0R) 

74  FORMAT  (3F20.0) 

55  IF(N0M)57,57,56 

56  READ  73,(XM(K),AMZ(K) ,K=1 ,NOM) 

57  IF(N0RM)59,5<3,58 

58  READ  74, (RXM(KA) ,RAMZ(KA),PSLP<KA) ,KA=1,N0RM) 

59  IF(N0I)61 ,61,60 

60  READ  72,(AE(  IA),BE(IA) ,EA(IA) ,EB(IA) , IA=1,N0I) 

61  IF(NCK)1 ,1,62 

62  READ  72, ( AG( IB) ,BG( IB), GA( IB) ,GB(IB ), IB=1,N0K) 
1  N0R1  =  NOR+1 

C1  =  0. 

C2=  0, 

N0RM2  =NORM  +  2 

N0RM3  =NORM  +  3 

NOQ  =  NOR  +NORM 

NCQ1  =  NOQ  +  1 

N0Q2  =  NOQ  +  2 
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N0Q3  =  NOG  +  3 

EPSX=10.**(-8) 

DO  2  K=1,NCR 

2  AA(1  ,K)  =  1. 

DO  3  K=  NORl,NOQ2 

3  AA( 1,K)  =  C. 
AAA=C. 

IF(NOF)U5,U5,35 
35  DO  U  N=l ,NCF 

k       AAA=  AAA  +  FY(N) 
45  IF(NOL)i*8,U8,U6 
U6  DO  5  1=1, NCL 

5  AAA  =  AAA  +  ( ( QA ( I ) +Q8 (  I  )  ) /2  .  )  *  (B(I)-A(I)) 
U8  AA( 1,NCQ5)=-AAA 

DO  6  K  =1,NCR 

6  AA(2,K)  =  C-XR<K) 
IF(NCRM)9,9,7 

7  DO   8  K=NORl,NOQ 

8  AA(2,K)=1. 

9  DO  1 C   K=NCQl,NOQ2 

10  AA(2,K)  =  0. 

CALL  MOMENT  (C,AAA) 
AA(2,N0G3)  =-AAA 
IF  (NORM) 17, 17,1  1 

11  DO  Id  J=3,N0RM2 
DO  13  K=1,NCR 
AA( J,K}=0. 

DO  12  IA=1,N0I 
AAA=EA( IA)*ENTON(RXM< J-2  ) ,AE(IA),XR(K),1)  -  EB(IA) *ENT0N(RXM (J-2) , 
1BE(  IA),XR(K),  1)  +((EB(IA)-EA(IA))/(BE( U)-AE(IA) > ) * { ENT1 N { RXM (J-2 ) 
2,AE(IA) ,XR(K) ,1  )-ENTlN(RXM(  J-2),  BE  {  IA),XR(K),1)) 

12  AA(J,K)  =  AA(J,K)  +  AAA 
IF  (NCK)  13,13,900 

900  DO  91C   IB=1,N0K 

IFtNORM-1 )902,9C1  ,9C2 
902  IF( J-NORM2J901 ,90U,90l 
9C1  X=RXM(J-2)  -EPSX 
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AAA  =  GA(IB)*  FTON(   X      , AG(  IB  ) , XR { K) , 0 )  -  GB(IB)*  FT0N(X 

1  ,BG( IB),XR(K) ,0)  +  ( (GB( IB)-GA(IB) )/(BG(IB)-AG( IB) ) )*  (  FT1N( 

2  X      ,AG(IE) ,XR(K),0)-  FT1N(   X      , BG (  IB  )  , XR( K ) , 0 ) ) 
GO  TO  910 

SOU  X=RXM(J-2)  +  EPSX 

AAA    =    GA(IB)*    FTON(       X  , AG ( IB ) , XR ( K ) , C )    -    GB ( I B ) »    FTON(X 

1  ,PG( IB) ,XR(K) ,0)    +    ( (GB( IB)-GA(IB))/(BG( IB)-AG(IB) J)*     (     FT1N( 

2  X  , AG(IB) ,XR(K) ,0)-    FTlNt        X  , BG ( I B ) , XR ( K ) , 0 ) ) 
910    AA{J,K)    =AA(J,K)     -AAA 

13  CONTINUE 
IF(N0C-N0R1 ) 1 55,  1  35,  135 

135    00    15    K=N0R1,NCG 
AA(J,K)=    0. 
DO    1U    IA=1,N0I 
AAA=    EAUA)«ENTON(RXM(  J-2  )  ,AE(  IA)  , RXM ( K-NOR ) , C )    -EB(IA)*    EKTCN 
1     (RXM (J-2),BE(IA)  ,RXM(K-NOR),0)     +     ( ( EB (  I  A )-EA ( I  A ) ) / ( BE { I  A )-A E (  I A  )  ) 
2)    *(ENTlN(RXM(J-2) ,AE( IA), RXM ( K-NOR) , 0 )    -    ENT 1 N ( RXK ( J-2 ) , BE (  I A ) , 

3  RXM(K-NCR) ,0) ) 

14  AA(J,K)     =    AA(J,K)     +    AAA 

15  CONTINUE 

155  AA( J.NCC1 )=1. 

AAU,N0C2)=    C. 
-       IF(J-N0RM2) 156, 157,156 

156  X=RXM(J-2)  -  EPSX 
CALL  SL0PE(X,AAA) 
GO    TO    16 

157  X=RXM(J-2)  +  EPSX 
CALL    SLCPE(X,AAA) 

16  AA(J,N0G3)     =-AAA     +    PSLP(J-2) 

17  DO    2C    J=N0RM3,N0G2 
DO    18    K=1,N0R 
AA(J,K)     =0. 

DO    131       IA=1,N0I 
AAA    =    EA(IA)«    DI TON(XR( J-N0RM2) ,AE(IA)  ,XR(K),  1)    - E B  ( I  A ) *C  IT CN 

1  (XR{ J-N0RM2) ,BE(  IA),XR(K) ,1 )     + ( { EB ( I  A )-EA ( I A )  ) /    { B E (  I A  )  -A E (  I A  )  )  ) 

2  *(OIT1N(XR(J-NORM2), AEt IA) ,XR(K),1 )    -    0  I  T 1 N ( XR ( J-NORf 2 ) , BE{  I A ) , 

3  XR(K),1) ) 
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181  AA(J*K)  =  AA(J,K)  +  AAA 

IF  (N0KJ18, 18,920 
920  DO  930   IB=1 ,NOK 

AAA    =GA(IB)»    ENT0NURU-N0RM2),  AG<  IB)  ,XR(K),  0)    -    GB(IB)«ENT0N 

1  (XR(J-N0RM2),BG(IB) ,XR( K) ,0)    + ( (GB ( IB )-G A( I B) ) / { BG(  IB)-AG( IB ) ) ) 

2  *(ENT1N(XR( J-NQRM2),AG(  IB)     ,XR(K),0)    -ENT1N ( XR ( J-N0RM2) , BG ( I  8) , 

3  XR(K),0) ) 

930    AA(J,K)    =AA(JtK)    -AAA 

18  CONTINUE 
IF(NOO-NORl) 195,185,185 

185    DO    19    K=N0R1 ,N0Q 
AA(J,K)    =    0. 
DO    19     IA=1,N0I 
AAA    =    EA(IA)«    DIT0N(XR( J-N0RM2)  ,  AE ( I  A) , RXM( K-NOR ) ,  0 )    -    EB(IA)* 

1  DITON(XR( J-N0RM2),BE( U ) ,RXM(K-NOR) ,0)    + ( ( EB( I  A  )-EA ( I  A) )     /     (BECIA 

2  )-    AE(IA)))     »    (DIT1N(XR( J-N0RM2) ,AE(IA) ,RXM(K-NOR) ,0)    -    DIT1N 

3  (XR{ J-N0RM2),BE(IA) ,RXM( K-NOR) ,0) ) 

19  AA(J,K)     =    AA(J,K)    +    AAA 
195    AA(J,N0Q1)    =XR(J-N0RM2) 

AA(J,N0Q2)     =    1. 

CALL    DEFLECT     ( XR( J-N0RM2 ) , AAA) 

20  AA(J,N0Q3)     =-AAA+    PDEFT ( J-N0RM2 ) 

"    CALL    GAUSS2(N0Q2, . 00000001, AA, XA,K1) 
GO    TO    (22,21),K1 

21  WRITE    OUTPUT    TAPE    4,77 

77       FORMAT    ( 16H    MATRIX    SINGULAR) 
STOP 

22  DO      23    K=1,N0R 

23  RY(K)  =  XA(K) 
IF(N0RM)26,26,24 

24  DO      25    K=N0R1 ,NOQ 

25  RAMZ(K-NOR)    =    XA(K) 

26  CI  =  XA(NOQl) 
C2  =  XAIN0Q2) 
X=0. 

NOPE    =N0P-1 
FNOP    =    NOP- 1 
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DO   27   J=1,NCPE 
CALL  LOAD  (X,W(J ) ) 
CALL  SHEAR  (X,V( J) ) 
CALL  MOMENT  U,AM(  J) ) 
XI  =X+EPSX 
CALL  SLCPE(X1,Y1( J)  ) 
Yl(J)=Yl(J)  +C1 

CALL  DEFLECT  (X, Y(J) ) 
Y( J)=Y< J)+C1«X+C2 
XX( J)  =  X 
27   X  =  C/FNOP  +X 
XX(NOP)=X 

X=X-3.*EPSX 
CALL  LCAD(X,W(NOP  )  ) 
CALL  SHEAR(X, V(NOP) ) 
CALL  MOMENT(X,AM(  NOP)  ) 
CALL  SLCnE(X,Yl (NOP)  ) 
Yl (NOP)=Yl (NOP5+C 1 
CALL  CEFLECT(X,Y( NOP) ) 
Y(NOP)=Y(NOP)+CHX  +C2 
WRITE  OUTPUT  TAPE  4,76 
76  FORMAT  (  3X  ,  1  HJ  ,7  X  ,4HL0AD,  10X,  5HSHEAR  ,  10X  ,  6HMCMENT  ,  8X  ,  5HSLCPE  , 
-  1  9X, 1CHDEFLECTI0N,8X,1HX) 

WRITE  OUTPUT  TAPE  4 , 78, ( J  ,  W ( J )  , V ( J ) , AM( J )  , Y 1  ( J )  , Y( J ) , XX { J )  ,  J  =  l  , 
1NOP) 
78  FORMAT(I4,6E15.8) 
IF  (NCR)282,282,280 

280  PRINT  281, (M,RY(M) ,M=1,NCR) 

281  FORMAT  (  3HRY(, 12 , 2H ) = ,E20. 8 ) 

282  IF(N0RM)288,288,283 

283  PRINT  284, ( KA, RAMZ ( KA ) , KA= 1 , NORM) 

284  FORMAT  ( 5HRAMZ ( , I  2, 2H )  =  , E20. 8 ) 
288  PRINT  286, CI 

286  FORMAT(5X,3HC1=,E20.8) 
PRINT  287, C2 

287  FORMAT(5X,3HC2=,E20.8) 

285  CONTINUE 
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STOP 
END 
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SUBROUTINE  LOAD(X,W) 

01  MEN  SI ON  A{25),B ( 25 ) f OA ( 25 ) , CB ( 25) , XFC( 100 ) ,FY ( 100 ) , XR ( 1 C ) , RY ( 1 C ) 
1,XM(5C),AMZ(50),RXM(10),RAMZ(  10),AE(25),BE(25),EA{25),EB<25) 
2,GA(25),GB(25),AG(25),BG(25) 

COMMON  A,BtQA,QB, XFOf FY,XR,RY,XM,AMZ»RXM, RAMZ,AE,BE,EA,ERiNOL, 
1  N0R,N0F,N0M,N0RMfNCI,N0K,GAtGB,AG,8G 

EPSX=10.*«(-8) 

w=o.  . 

1  DO  5  I*1,N0L 

2  IF(X  +  EPSX-B(I) )3,  3, 5 

3  IF(X+EPSX-A(I))5»4,1j 

k    WW=QA(I)+((QB( I )-QA{ I))»(X-A(I) ))/(B( I)-A(I ) ) 

W=W+WU 
5  CONTINUE 
RETURN 
END 


5^ 


SUBROUTINE  SHEAR(X,V) 

DIMENSION  A(25),B(25) ,QA ( 25 ) , QB ( 25) ,XFC(  IOC ) , FY ( 100 ) , XR ( 10) ,RY(10) 
1,XM(5C) ,AMZ(50),RXM{1C),RAMZ(10) ,AE(25)fBE(25),EA(25),EB(25) 
2,GA{25),GB(25),AG(25),BG(25) 

COMMON  A,B,QA,QB,  XFO, FY, XR, RY , XM, AMZ, RXM, RAMZ t AE, B E» E A, EB, NOL , 
1  N0R,NOF,NOM,N0RM,NCI,N0K,GA,GB,AG,BG 

EPSX=10.»»(-8) 

V=0. 

IF(N0L)7,7,1 

1  DO  6  1=1, NOL 
IF(X-Ad)  )6,2,2 

2  IF(X+EPSX-B(I ) )U,  4,3 

3  VV=  QA(I )»(X-A(I)  )  +  (  (QBU  )-QA(I ) ) • ( (X-A t I ) )**2) ) /<  2 .«( B( I )-A ( I ) ) ) 
1-QBU  )*(X-B(  I)  )  - ( (QB(I)-QA(IJ ) *( (X-B I  I ) ) »»2) )/ (2. » <B ( I >-A ( I ) )  ) 

GO  TO  5 
U    VV  =+QMI)*(X-A(  I)  )♦(  (QB<  I  )-QA(I  )•)•(  (X-A(  I))»*2))/(2.*(B(I)-A(I)J) 

5  V=V+VV 

6  CONTINUE 

7  IF(N0F)10,10,75 
75  DO  9  N=1,N0F 

IF(X+EPSX-XFO(N))  9,9,8 

8  VV=FY(N) 
V=V+VV 

9  CONTINUE 

10  IF(N0R)13,13,105 
105  DO  12  M=l ,NOR 

IF(X+EPSX-XR(M) )12,12,11 

11  VV=  RY(M) 
V=V+VV 

12  CONTINUE 

13  CONTINUE 

RETURN 
END 
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SUBROUTINE  FOMENT ( X, AM ) 

DIMENSION  A(25),B(25),QA(25),QB(25) ,XFO(100) ,FY(100),XR(1C),RY(10) 
1,XM(5C) ,AMZ(50) ,*XM(10),RAMZ(10),AE(25),BE(25) ,EA(25),EB(25) 
2,GA(2  5),GB(25),AG (25) ,BG(25) 

COMMON  A,BtQA,QB, XFO , FY ,XRt RY, XM, AMZ, RXM ,RAMZ , AE ,B E, EA, EB, NCL , 
1  NORt  NOFt  NOMfNORMffNOItNOKf  GAf'GBt  AGtBG 

EPSX=10.««(-8) 

AM=0. 

IF(NOL)8,8,2 

2  DO  7   1=1 ,NOL 
IF(X-M  I)  )7,3,3 

3  IF(X  +  EPSX-B(I) )5,  5,4 

4  AMM=( +QA(I)/2.)*( (X-A(I) )««2)+( (QB(I)-QAU) )*( (X-A(I) )*»3) J/(6.*( 
1 B( I )-A(  I)  ) )-(OB(I  )/2.)*( (X-B( I) )»«2>-  (QB(  I  )-CA{  I )  )« ( ( X-B ( I )  )**3  ) 
2/(6.»(B(I )-A( I))) 

GO  TO  6 

5  AMM=(+QA( .-!  )/2.)»(  (X-A(  I)  )**2)+(  ( QB ( I ) -QA ( I ) ) * ( (X-A ( I ) ) *»3 ) ) / ( 6. * 
1 (B(  I)-A(I  )  >  ) 

6  AM=AM+AMM 

7  CONTINUE 

8  IF(NOF)12,12,9 

9  DO  11  N=l ,NOF 

-  IF(X+EPSX-XFO(N))  11,11,10 

10  AMM=  FY(N)*(X-XFO (N) ) 
AM=AM4AMM 

11  CONTINUE 

12  IF(NOR) 16, 16,13 

13  DO  15  M=l ,NOR 
IF(X*EPSX-XR(M) )1 5,15,14 

14  AMM=  RY(M)*(X-XR(  M)  ) 
AM=AM+AMM 

15  CONTINUE 

16  IF(N0N)28,28,  17 

17  DO  19  K=l ,NOM 
IF(X+EPSX-XM(K))1 9,19,18 

18  AMM=ANZ(K) 
AM=AM+AMM 
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19  CONTINUE 

28  IF<NORM  32,32,29 

29  CO  31  KA=1,NORM 
IF<X+EPSX-RXM(KA) )31,31,30 

30  AMM=R/WZ(KA) 
AM=AM+AMM 

31  CONTINUE 

32  CONTINUE 
RETURN 
ENO 
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SUBROUTINE    SLOPEl X,Yl  ) 

DIMENSION    A(25),B (25) ,QA{25) ,QB(25),XFC(  100 ) , FY ( 100 ) , XR ( 10) , RY {  10) 
ttXM(5C),AM2(50),RXM( 1 0 ) , RAMZ( 10 ) , AE (2  5 ) ,BE (25 ) ,E A( 25) ,EB(2  5) 
2fGA(25),GB(25),AS(25),BG(25) 

COMMON    A,B,QA,GB, XFO, FY, XR,RY, XM, AMZ, RXM ,RAMZ , AE,B E, EA, EB, NOL  , 

1  NOR, NOF,NOM, NORM tNOI ,NOK ,GA , GB , AG, BG 
Y1  =  0. 

DO    17       IA=1,N0I 
IF(X-AE( IA)  )  17, 17  ,2 

2  IF(N0L)5,5,3 

3  DO      4       1=1, NOL 

20  EYY    =    EA(  IA)*( (QA ( I ) /2. ) *ENTON( X, AE ( I  A) ,A( I) ,2)     +  (  ( GB ( I )-QA( I )  ) / 

1 (6.*(B( I )-A(I ) ) )) *(ENTON(X,AE( I  A) ,A(I ) , 3 J-ENTCN ( X, AE( I A ) ,B ( I } , 3 j ) 

2  -(Q8(I)/2. )»ENTON(X, AE( IA),B(I),2))       +  (  ( EB ( I  A )-EA ( I  A) ) / { BE (  I A  )~ 

3  AE(IA) ) )*( (QA(I) /2.)»ENT1N(X,AE( IA),A(I ),2)    +  ( (GB { I ) -CA (  I  )  )     / 
4(6.*(E(I)-A( I ) ) ))*(ENT1N(X,AE(IA)  , A ( I ) ,3 )~ENT IN { X,  AE ( I fi ) t 8 ( I  ) , 3)  ) 
5    -(QB( I )/2. )*ENT1N(X,AE( IA),B(I),2)) 

Y1=Y1+EYY 

21  EYY=    -EB(  IA)»( ( QA ( I ) /2. ) *ENTON( X, BE ( I  A) ,A(I ) ,2)    + (  ( GB ( I )-GA(  I  )  )  / 
1(6.»(B( I) -A (I ) ) )) «(ENTON(X,BE(IA) , A( I ) , 3 )-ENTON (X, BE (  I  A  ) , 6 ( I  )  , 3 ) ) 

2  -(QB(I)/2.)»ENT0N(X,BE( IA),B(I),2))       -(  (EB(IA)-EA(  I A  )  )  / ( BE (  I A  )- 

3  AE(IA) ) )*( (QA(I)/2.)»ENT1N(X,BE( IA), A(I ),2)     + ( ( GB ( I )-QA ( I ) )     / 
**(6.*(E(  I)>-A(I  )  )  ))«(ENT1N(X,BE(IA)  ,A(I)  ,3  )-ENT  1  N(  X,  BE  (  IA),B(I),3)) 
5    -(QB(I )/2. )*ENT1N(X,BE( IA),B(I),2)) 

Y1=Y1+EYY 
U    CONTINUE 

5  IF(N0F)8,8,6 

6  DO       7    N=1,N0F 

22  EYY=  EA( IA)«FY(N) «ENTON(X,AE( IA ) , XFO( N) , 1 )  -EB ( I  A) *FY (N ) *ENT0N ( X  , 
1BE( IA),XFO(N), 1)  +  ( (EB( IA)-EA(IA) )/(BE( IA)-AE(IA) ) )«FY(N)« 
2(ENT1N(X,AE(IA),XF0(N),  1)  -  ENT1 N ( X,BE ( I  A  )  ,XFO (N ) ,  1 ) ) 

Y1=Y1+EYY 

7  CONTINUE 

8  IF(NOR)  11,  1  1,9 

9  DO  10  M=1,N0R 

23  EYY=  EA(IA)«RY(M) *ENTON( X , AE ( I  A) , XR ( M ) , 1 )-EB( IA)*RY(M)»  ENTCN(X, 
1BE(IA),XR(M),  1)  + ( (EB( IA)-EA(IA) ) / ( BE ( I A)-AE (  IA) ))*RY(M)  « 

58 


2(ENT1N(X,AE(  IA),XR(M),1)-    ENT1 N( X, BE(  I A )  ,  XR (M ) , 1 ) ) 
Y1=Y1+EYY 

10  C0NTTNUE 

11  IF(NOM     14,14,12 

12  00    13      K«lfNOM 

2U  EYY=  EAUA)«AMZ(K)«ENTON(X,AE(IA),XM(K),0)  -EB( I  A) »AMZ < K) »ENTON (X , 
IBE(IA)  ,XP(fc),0)  +  <  (EB(IA)-EA(IA)  )/(BE<  IA)-rAE(IA))  )•  AMZ(K)  • 
2  (ENT1N(X,AE( I A ) , XM ( K) ,0)-  ENT1 Nl X*BE ( I A ) , XM < K) ,0)  ) 
Y1=Y1+EYY 

13  CONTINUE 

14  IF(NORM)  17,17,15 

15  DO  16   KA=  1,N0RM 

25  EYY  =  EA( I A )»RAMZ(K A)«ENTON(X, AE( I  A ) , RXM (KA) ,0)-EB ( I  A )*RAMZ ( KA ) « 

1  ENTON(X,BE(IA),RXM(KA),0)  +  < ( EB < I A) -E A( I  A ) ) / ( BE( I A)-AE(  I A  )  )  )  » 
2RAMZ(KA)»(ENT1N(X,AE( I A ) , RXM ( KA ) ,0 )-  ENT IN (X,BE ( I A ) ,RXM<KA ) ,0 )  ) 
Y1=Y1+EYY 

16  CONTINUE 

17  CONTINUE 
IF(NCK)595, 595,500 

500    DO    59C       IB=1,N0K 

IF    (X    -AG(IE) )590,59C,509 

509  IF    (NCL)530,530,510 

510  00    52C       1=1, NOL 

26  EYY    =    GA(  IB)»(QA(  I  )»FTON  (  X,  AG(  I  B)  ,A(I  ),1  )    ♦    ((QB(I')    -QA(I)     )    / 

1  (2.«(B( I)-A( I ) )) )»(FTON(X,AG(IB) ,A(I ) , 2 )-FTON( X ,AG( I B ) , B ( I ) ,2)) 

2  -Q8( I)»FTON(X, AG ( IB) ,B(I),1) )    ♦    ( ( GB ( IB )-GA ( I B) )/ ( BG( IB )-AG ( I B  )  ) ) 

3  «(QA(I )«FT1N(X,AG(IB),A(  I),l  )    +    ( (QB ( I )-QA( I ) )/( 2. • ( B ( I )-A(  I )  )  )  ) 
4»    (FT1N(X,AG(IB), A(I),2)    -FT 1N( X, AG( I B) , B( I ) , 2 )  )    -    QBU)»    FT1N 

5    (X,AG(  IB),B(I) ,2) ) 
Y1    =Y1    -EYY 

27  EYY    =-GB(  IB)*(QA(  I  ) «FTON( X, BG( IB ) , A( I  ),  1)    ♦    (<QB(I)    -CA(D)    / 

1     <2.»(B( I)-A( I ) )) )*(    FTON(X,BG(IB) ,A{ I ) , 2 )-FTON (X, BG ( I B ) , B ( I ),2) ) 
2-QB{ I )•    FTON(X,BG( IB),B(l),l))-    ( (GB( I B )-GA ( IB) )/(BG( IB)-AG( IB)) ) 
3*(QA(I)»    FT1N(X,BG(IB),A(I),1)+    (  (  QB(  I  ) -QA ( I  ) )/ ( 2.  • (B ( I )-A < I  )  )  )  ) 
4»(     FT1N(X,BG(IB),  A(I),2)-    FT  1  N(  X,  BG(  I  B)  ,  B(  I  )  ,2) )    -    QB(I)«      FT1N 
5    (X,BG( IB),B(I),2)) 
Yl    =Y1    -EYY 
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520  CONTINUE 

530  IF(NOF)56C,560.540 

5U0  DO  55C   N=l,NOF 

28  EYY*GA< IB)*FY(N)»  FTON (X , AG( I B ) , XFO( N ) ,0 )  -GB ( I B )* FY ( N )«  FT0N(X, 
1BG(  I8),XF0(N)  ,0)  +  (  (GB,(  IB)-GA(IB)  )  /  (  BG(  IB  ) -AG  ( I  B  )  )  )  *  FY(N)* 
2(  FT1N(X, AG( IB).XFO(N) ,0)  -   FT1 N( X,BG( I B ) ,XFC (N ) , C ) ) 

Yl  =Yl  -EYY 
550  CONTINUE 
560  IF(N0R)59C,590,570 
570  DO  58 C   !!-l,NOR 

29  EYY  »GA(IB)*RY(M)«  FTON ( X , AG( I B ) , XR(M ) ,0  )  -GB  ( I  B  )*  RY  (  M  >  *  FTONtX, 
1BG(IB),XR(M), 0)  +  { (GB(IB)-GA(IB) )  /  ( BG ( I B J-AG ( IB ) ) )  *  RY(M)  * 
2(  FT1N{X» AG(IB) ,XR(M), C)  -   FT1  N ( X, BG ( I B ) , XR (M) , 0)  ) 

Yl  =  Yl  -EYY 
580  CONTINUE 
590  CONTINUE 
595  CONTINUE 

RETURN 

END 
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•  SUBROUTINE  DEFLE: T ( X, Y ) 

01  MENS  I ON  A (25) ,B (25) ,QA ( 25) , CB ( 25 ) ,XFO ( 100) , FY( 100 ) , XR( 10 )f RY(10) 
1,XM(5C),AMZ(50),RXM(10),RAMZ( 10 ) , AE(25) , BE (25 ) , EA( 25 ) , EB( 25 ) 
2,GA(25),GB(25),AG(25),BG(25) 

COMMON  A,B,QA,QB,  XFO, FY, XR,RY ,XM, AMZ, RXM, RAMZ , AEtBEt EA ,EB,NOL, 

1  NOR,NOF,NOM,NORM,NOI,NOK,GA,GB,AG,BG 
Y=0. 

DO  17   IA=1,N0I 
IF(X-AE(IA))17, 17,2 

2  IF(N0L)5,5,3 

3  00   4   1=1, NCL 

20  EYY  =  EA( IA)*( ( QA ( I )/2. ) *DITON( X, AE( I  A) , A ( I ) , 2)  ♦  (  ( QB ( I )-QA (  I )  )/ 

1 (6.«(E(I)-A(I) ) )) •(DITON(X,AE(IA) ,A(I ) , 3 )-DITON(  X,  AE ( I  A) ,B( I),3)) 

2  -(GB(I)/2.)«DIT0N(X,AE(IA),B(I),2))   ♦  ( ( EB ( I  A )-EA ( I A ) )/ ( 8E( I A)- 

3  AE( IA) )  )*( (QA(  I)/2.)»DIT1N(X,AE(  IA), A( I),2)  +  (  (QB  (  I  )-QA  (  I  )  )  / 
U(6.*(B(I)-A(I)) ))»(DIT1N(X,AE(IA)  ,A ( I ) ,3 )-DIT 1N( X,  AE(IA)  ,B(I)f3)) 
5-(QB(I)/2.)«DIT1N(X,AE(IA),B(I) ,2) ) 

Y=Y+EVY 

21  EYY=  -EB( IA)»((QA(I)/2.)«DIT0N(X,BE(IA),A(I),2)  ♦  (  ( QB ( I )-QA(  I )  )/ 

1 (6.»(B( I) -A (I ) ) ))«(DIT0N(X,BE(IA),A(I ) ,3) -DITON( X,  BE ( I A ) , B ( I ),3) ) 
2-(QB(I)/2.)»DIT0N(X,BE(IA),B(I),2))   -( ( EB ( I  A) -EA ( I A ))/ ( BE ( I A ) - 
3AE(  IA)n*((QA(I)/2.)«DITlN(X,BE(.IA),A(I),2)  -M(QB(  I)-QA(I))  / 
4(6.»(B(I)-A(I)))) *(0IT1N(X,BE(IA),A(I),3)-DIT1N(X,  BE (I  A) ,B(I),3) ) 
5-(QB(I)/2.)«DIT1N(X,BE(IA),BU),2)) 
Y=Y+EYY 
U  CONTINUE 

5  IF(N0F)8,8,6 

6  DO   7  N=1,N0F 

22  EYY=  EA( IA)*FY(N) «D I TON( X , AE ( I A ) , XFO( N  ) , 1)  -EB ( I  A) «FY( N) »D ITON (X, 
1  BE(IA),XFO(N), 1)  ♦  ((EB(IA)-EA(IA))/(BE( I  A)-AE ( IA ) ) ) »FY (N )• 
2(DIT1N(X,AE(IA),XF0(N),1)  -  DIT IN (X ,BE( I A ) , XFO( N ) , 1 ) ) 

Y=Y+EYY 

7  CONTINUE 

8  IF(N0R)11,U,9 

9  DO  10  M=l ,NOR 

23  EYY=  EA(IA)*RY(M) »DITON( X, AE ( I  A ) ,XR ( M  ) , 1 )-EB (  IA ) *R Y( M )•  DITON(X, 
1BE( IA),XR(M),1)  +((EB(IA)-EA( IA) ) /( BE ( I A )-AE ( I  A) ))»RY(M)  • 
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2{DIT1N(X,AE(IA),XR(M),1)-  DITlN(X,BE(IA),XR(M),l)) 
Y=Y+EYY 

10  CONTINUE 

11  IF(NON) 14,14,  12 

12  DO  13   K=1,N0M 

24  EYY=EA( IAJ#AMZ(K) «DI TON ( X , AE ( I A ) , XM( K ) , 0 )  -EB( I  A ) • AM2 ( K ) »D I TON ( X , 
1BE( IA),XM(K),0)  +  ((EB(IA)-EA(IA)  )/(BEUA)-AE(IA))  )*AMZ(K)* 

2  (DIT1N(X,AE( IA),XM<K),0)-  DIT  IN < X,BE ( I A ) , XM{ K) ,0)  ) 
Y=Y+EYY 

13  CONTINUE 

14  IF(NORM)  17,17,15 

15  DO  16   KA=  1,N0RM 

25  EYY=EA(IA)»RAMZ(KA)*OITON(X,AE(IA) , RXM(KA) ,0 ) -EB ( I  A ) «RAMZ ( KA  )» 

1  DITON(X,BE(IA),RXM(KA),0)  +  ( ( EB( I  A) -EA ( I A ) ) / { BE( I  A) -AE (  I A )  )  )  « 
2RAMZ(KA)*{DIT1N(X ,AE(IA) ,RXM(KA),0)-  DIT 1 N( X, BE ( IA ) ,RXN( KA  ) , C  )  ) 
Y=Y+EYY 

16  CONTINUE 

17  CONTINUE 
IF(NOK)595,595,500 

500  DO  59C   IB=1,N0K 

IF  (X  -AG( IB)  )590,590,509 
509  IF  (NCL)530,53Ct5  10 
51CT  DO  52C   .I^l.NOL 

26  EYY=GA(IB)*(OA{ I) »  ENTON ( X , AG { IB) , A{ I ) , 1 )  +  ((QB(I)  -GA(I)  )  / 

1  (2.*(B( I)-A( I) ))  )*( ENTON {X, AG ( IB) ,A( I) ,2  )-ENTON (X f AG (  IB),B(  I  ),2)  ) 
2-QB(I)*ENT0N(X,AG(IB) ,B( I ),1) )  +  ( ( GB ( IB)-GA( IB) )/ (BG( IB )-AG ( I B) ) ) 
3*(QA( I)*ENT1N(X,AG(  IB),A( I  )  ,1 )  +  { (QB( I )-CA{ I ) ) / (2 . * ( B ( I ) -A {  I )  )  )  ) 
4*(ENT1N(X,AG( IB),  A(I),2)-ENT1N(X,AG(IB) ,B( I),2) )  -  G8(I)*ENT1N 
5  (X,AG(IB),B( I), 2  )  ) 
Y=Y-EYY 

27  EYY=-GB( IB)*(QA(I  ) »ENT0N ( X ,BG { IB) ,A( I  ), 1  )  + ( ( QB { I ) -QA (  I )  )  / 

1 (2.»(B( I )-A { I ) ) ))*(ENT0N(X,BG( IB),A(I ) ,2 ) -ENTON ( X, BG ( I B )  ,B(  I), 2)  ) 
2-QBU  )« ENTON  (X,BG(  IB  )  ,  B(  I  )  ,  1  )  )-  (  (GB{  IB)-GA(  I  B  )  )  /(  BG  (  I  B  )-AC  (  IB  )  )  ) 
3*(QA( I)«ENT1N(X,BG( IB),A( I ), 1)+  ( (QB( I )-QA ( I ) ) /( 2. •( B ( I )-A ( I ) ) ) ) 
4»(ENT1N(X,BG( IB), A(I ) , 21-ENT1 N( X, BG( I B) , B(  I), 2)  )  -  GB(I)«  ENT1N 
5  (X, BG(IB) ,B(  I), 2  )  ) 
Y=Y-EYY  ! 
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520  CONTINUE 

530  IF(N0F>560,560,540 

5U0  DO  550   N=1,N0F 

28  EYY*GMIB)»FY(N)*ENTON(X,AG(  IB)  ,XFO(N)  ,0)  -GB  (  IB  )*F  Y(  N  )  »ENTCN  (  X, 
1BG( IB)fXFO(N) f0)  +  UGB(  IB)-GA(IB)  )  /  (BG ( IB ) -AG ( I  B ) ) )  •  FY(N)* 
2(ENT1N(X,AG( IB),XFO(N) ,0)  -  ENT 1N( X,BG ( I B ) , XFO ( N)  ,  C) ) 

Y=Y-EYY 
550  CONTINUE 
560  IF(N0P)590, 590,570 
570  DO  580   M* 1, NOR 

29  EYY  =GA{ IB)*RY(M) *ENTCN( X,AG(IB) fXR(M),0)  -GB ( IB )*RY(M)*ENT0N( X, 
1BG( IB)tXR(M)tO)  +  ( (GB( IB)-GA(IB) )  /  ( BG ( I B)-AG ( IB ) ) )  «  RY(M)  * 
2(ENT1N(X,AG(IB),XR(M) ,0)  -  ENT 1N( X,BG ( IB ) , XR ( N ) ,0)  ) 

Y=Y-EYY 
580  CONTINUE 
590  CONTINUE 
595  CONTINUE 

RETURN 

END 
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FUNCTICN  FTON(XfA,B,N) 
FTON  =  0. 
IF(X-A)3,3, 1 

1  IF(X-e)3,3,2 

2  FTON=(X-B)*«N 

3  RETURN 
END 


FUNCTION  FT1N  (X,A,B,N) 

FT1N  =  0. 

IF  (X-A)3,3,l 

1  IF  <X-B>3,3,2 

2  FT1N  =  (X-B)**{N+1)  +  (B-A)»( (X-B)**N) 

3  RETURN 
END 
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FUNCTION   ENTCN  (X,A,B,N) 
FN=N 

ENTON  =  0. 
IF(X-A)5,5,1 

1  IF(X-e)5,5,2 

2  IF(A-B)3,3,U 

3  ENTON  =  (X-B)«»( N+1)/(FN+1.) 
GO  TC  5 

4  ENTON  *  (X-B)»«(N+1 )/(FN+l.)  -  ( A-B) «•( N+ 1 )/( FN+1 . ) 

5  RETURN 
END 


FUNCTICN   ENT1N( X,A, B,N) 
FN=N 

ENT1N  -    0. 
IF(X-A)5»5,1 

1  IF(X-e)5,5,2 

2  IF(A-B)3,3,4 

3  ENT1N  =  (X-8)«*(N+2)/(FN  +  2.  )  +   ( B-A )»( X-B )*•( N+l )/( FN+1 . ) 
GO  TQ  5 

U   ENT1N  '    (X-B)««(N+2)/(FN+2.)  +  { B-A) «<X-8) •* (N+l )/( FN+1 . ) 

1  -(A-B)«*(N+2)/(FN+2.)  -(( B-A )» (A-B) ♦«( N+ 1  ))/( FN  +  1  .  ) 
5   RETURN 
END 
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FUNCTION  DITONJX,  A,B,N) 
FN=N 

DITON  =0. 
IF(X-A)5,5, 1 

1  IF(X-e)5,5,2 

2  IF(A-E)3,3,U 

3  OITON  =  (X-B)«»(N+2)/((FN+l ,)#(FN+2.) ) 
GO  TO  5 

4  DITON  =  (X-B)*»(N+2>/( (FN+l.)«(FN+2.) )  - (X* ( A-9 ) «* ( N+ 1 ) ) / ( FN+1 . ) 
l-(A-B)«MN+2)/((FN+l.)«(FN+2.))   +  ( A* ( A-B ) »» ( N+1 )  > / ( FN+ 1 .  ) 

5  RETURN 
END 


FUNCTION   DIT1N(X,A,B,N) 

FN=N 

DITlN^O. 

IF(X-A)5,5,1 

1  IF(X-e)5,5,2 

2  IF( A-E)3,3,U 

3DIT1N  =tX-B)*»(N+3)/( ( FN+2. )*( FN+3 . ) )  +  ( ( B-A )» { X-B) *« (N+2 ) ) / 
1((FN+l.)»(FN+2.)) 
GO  TO  5 

4  DIT1N  =  (X-B)»»(N+3>/( ( FN+2. )*( FN+3. ) )  +  ( ( B-A ) » (X-B ) ♦» ( N+2 ) ) / 
l((FN+l.)*(FN+2. ))-(X«(A-B)»»(N+2) )/(FN+2. )-( X* (B-A ) » ( A-B ) ** ( N+ 1 ) )/ 

2  (FN+1.)  -  (A-B)»»(N+3)/( ( FN+2. )♦( FN+3. ) )  -  ( B-A )» ( A-B )»* ( N+2 ) / 

3  ( (FN+1.)»(FN+2.) )  ♦  (A»(A-B)»*(N+2))/(FN+2.)  +  ( A* ( B-A) * ( A-B ) ** 

4  (N+1))/(FN+1.) 

5  RETURN 
END 
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SUBROUTINE      G AUSS2 ( N, EP, A, X, KER) 
(DIMENSION    A<22,23)     ,X(22) 
NPM=N+1 

10  DO    31*    L=1,N 
KP«0 

z-o.c 

DO    12    K=L,N 
IF(Z-ABSF(A(K,L>  ))11,12,12 

11  Z=ABSF(A(K,L)) 
KP=K 

12  CONTINUE 
IFU-KPJ13,  20,20 

13  DO  1*4  J=L,NPM 
Z=A(L,J) 

A(L, J)=A(KP,J) 
1U   A(KP,J)=Z 
20   IF(AeSF(A(L,L) )-EP  )50,5C, 30 

30  IF(L-N)31,40,40   ' 

31  LP1=L+1 

DO  314  K*LP1,N 
IF(A(K,L)>32,34, 32 

32  RATIO=A(K,U/A(L,L) 
DO  33  J=LP1,NPM 

33  A(K,J)=A(K, J)-RATIO»A(L, J) 
3U  CONTINUE 

40  DO  43  1=1, N 
II=N+1-I 

DO  43  J=l,l 

JPN=J+N 

S*0.0 

IF(II-N)41,43,43 

41  IIP1-II+1 

DO  42  K=IIP1,N 

42  S=S+A(II,K)«X(K) 

43  X(II)   =  (AUI,JPN)-S)/A(II,II  ) 
KER=1 

75    RETURN 
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50   KER=2 
RETURN  . 
END 
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